# Wilcoxon matched-pairs signed-ranks Test
# By Tecumseh Fitch. Version 0.4
# This code has been tested against both test data sets in Siegel and matches his results exactly
# (though note a mistake in his Table 5.6 - last data value should be 62)
# Calculation of p values uses code tested in "NormalDistribution.py" and borrowed from JSci

from math import sqrt, exp

def erfc(x):
    """complementary error function=2/sqrt(pi) \int_x^{\infty} dt e^{-t^2}. Code from numerical recipes."""
    z=abs(x)
    t=1.0/(1.0+0.5*z)
    ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))))
    if x>=0.0:
        return ans
    else:
        return 2.0-ans
    
def CumulativeNormalDistribution(z):
    """cumulative normal distribution probability, with mean 0 and variance 1"""
    return 1.0-0.5*erfc(z/sqrt(2.0))

def sign(x):
    if x < 0:
        return -1
    else:
        return 1

def tiedRanks(val):
# Takes a sorted list of values, and outputs their ranks dealing with ties
# A string of tie values are replaced with their average rank, so 1, 2, becomes 1.5, 1.5
# Give ranks starting with 0! Note that Siegel starts ranks at 1
        
    rank= [x for x in range(len(val))] # generate the preliminary list of ranks, without tie adjustment

    # Now adjust for ties: A surprisingly tricky problem!
    
    first = total = 0 # the trick for this algorithm is that the first duplicate, and the initial sum
                        # are constantly being set, and then thrown away if there is no duplicate.
    
    N = len(val)
    run = False
    for i in range(0, N-1):
        if val[i] != val[i+1]: # if the current two numbers are different, we are either at the end of a string of ties, or home free
            if run == False: # no ties yet!
                first = total = i+1 # store the first rank, and add to the sum
            else: # this _was_ the end of a run of duplicates: replace them all with the tie value
                run = False
                last = i
                average = total/(last-first+1.0)
                #print "replacing string of duplicates from", first, "to", last, " with average of ", average
                for j in range (first, last+1):
                    rank[j] = average
                first = total = i+1
        else: # we have a duplicate!
            run = True
            total += i+1 # add the duplicate's rank to the sum (first is already there)
    # end of the loop: but what if the very last value was a duplicate?
    
    if run == True: # the final value was a duplicate: current "first" value and total are OK
        last = i+1
        average = total/(last-first+1.0)
        #print "replacing string of duplicates from", first, "to", last, " with average of ", average
        for j in range (first, last+1):
            rank[j] = average

    return rank

def Wilcoxon(x, y):
# Implements the Wilcoxon matched-pairs signed-ranks test, each data point in x paired with one in y
# Basically follows Siegel 1965
# Note that the output statistic T (also called W by some) retains its sign, Siegel gives it unsigned


    returnArray=[-1,-1,-1]

    if len(x) != len(y):
        print "Error: data lengths must match"
        return returnArray
    if len(x)==0:
        print "Error: empty data set"
        return returnArray
    
    print "x:", x, "y:", y
    diffs = []
    for i, xval in enumerate(x):
        diffs.append(xval - y[i])

    print "diffs:", diffs

    # First delete differences of zero:
    diffs[:] = [x for x in diffs if x !=0]
    print "diffs omitting zeros:", diffs

    # Sort 'diffs' by their _absolute_ value, but retaining the sign:
    diffs.sort(lambda x,y:cmp(abs(x),abs(y)))
    #diffs.sort
    print "sorted:", diffs

    #generate a new list with the absolute values of the differences:
    absdiffs = [abs(x) for x in diffs]
    ranks = tiedRanks(absdiffs)
    ranks = [x+1 for x in ranks] # add 1 to every rank to adjust for the fact that ranks in tiedRanks start at 0, not 1

    # Re-assign the original signs to the ranks
    signedRanks = []
    for i, x in enumerate(ranks):
        signedRanks.append(x*sign(diffs[i]))
    print "signed ranks:", signedRanks
    
    negRanks = [x for x in signedRanks if x < 0]
    negT = sum(negRanks)
    posRanks = [x for x in signedRanks if x > 0]
    posT = sum(posRanks)

    sumrankssigned = sum(signedRanks)
    sumranksunsigned = sum(ranks)
    sumsqranks = sum([x*x for x in ranks])
    denom = sqrt(sumsqranks)
    
    print "sum (all signed ranks)=" , sumrankssigned, "sum/2=", sum(signedRanks)/2.0
    print "sum (unsigned ranks)=", sumranksunsigned
    print "T (Unk version)=", sumrankssigned/denom
    print "negative T:", negT, "positive T:", posT

    if abs(negT) > abs(posT):
        T = posT
    else:
        T = abs(negT)
        
    # finally calculate the z score (formula 5.5, pg 79 in Siegel 1965 "Nonparametric statistics")
    N = len(diffs)
    m = (N*(N+1))/4.0
    sd = sqrt((N*(N+1)*(2*N+1))/24.0)
    
    #what to do if T-m is 0.0???
    # it means that there is no difference in the sample
    
    
    if (T-m) != 0.0:
        z = (T - m)/sd
        zcont = (T - m - 0.5)/sd
        print "N =", N
        print "Wilcoxon T Statistic=", T
        print "z score =", z
        print "z score with continuity correction =", zcont #(formula 2.2, p 17 Gibbons (1993) "Nonparametric statistics")
        print "p (uncorrected, one-tailed)=", CumulativeNormalDistribution(z)
        print "p (uncorrected)=", CumulativeNormalDistribution(z)*2.0
        print "p (corrected)=", CumulativeNormalDistribution(zcont)*2.0   
        if N < 25:
            print "Warning! Z-score calculation may be inappropriate for N less than 20 or 25!"
        returnArray[0] = CumulativeNormalDistribution(zcont)*2.0
        returnArray[1] = z
    else: 
        returnArray[0]=1.0
        returnArray[1]=0.0
    returnArray[2] = T
    
    return returnArray

# To run this do something like:
##x = [1, 2, 3, 4, 5, 6, 7, 8]
##y = [1, 3, 2, 6, 8, 9, 7, 10]
##Wilcoxon(x, y)
